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SUMMARY 

This  document  reports  on  the  work  undertaken  in  the  first  9  months  of  the  project. 


BACKGROUND 

The  aim  of  this  project  is  to  undertake  a  feasibility  study  into  the  potential  utility  of 
integrating  high  resolution  two  dimensional  finite  element  flow  models  and 
Geographical  Information  Systems  technology. 

The  initial  phase  of  this  research  concerns  the  construction  of  an  operational  high 
resolution  flow  model  for  a  60km  reach  of  the  Missouri  River  between  Gavins  Point 
Dam  and  Maskell  gauging  station.  Specifically,  the  contract  seeks  to  produce  a  CRREL 
report  and  to  assess  data  needs  for  2D  FE  models  for  river  flow  inundation. 

This  report  contains  a  brief  review  of  progress  on  this  work  unit  during  months  6-9  of 
the  research  contact. 


PROGRESS 

This  reporting  period  saw  the  production  of  the  draft  CRREL  report  on  the  contract 
—  working  title  ‘’'’Towards  an  evaluation  ofl&  2D  models  for  floodplain  inundation 
mapping'. 


A  copy  of  the  draft  report  is  attached. 
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1.  INTRODUCTION 

In  this  report  we  address  the  numerical  simulation  of  reach  scale  free  surface.  For  such  flows 
construction  of  accurate  prediction  models  remains  a  fundamental  research  problem  in 
hydrology  and  hydraulics.  Such  studies  are  not  only  of  importance  to  engineers  and  planners 
engaged  in  the  water  resources  field  but  are  also  an  important  pre-requisite  for  modelling 
investigations  of  sediment  transport  (James,  1985;  Pizzuto,  1987;  Marriot,  1992;  Falconer 
and  Chen,  1996;  Walling  et  al.,  1996),  contaminant  transport  (Marron,  1989;  Burt  et  al., 
1993),  floodplain  ecology  (Jolly,  1996;  Richards  et  al.,  1996)  and  catchment  hydrology 
(Beven,  1985;  Abbott  et  al.,  1986;  Beven  et  al.,  1987)  which  either  require  hydraulic  data  as 
an  input  or  incorporate  a  submodel  for  hydraulic  processes.  For  all  these  tasks  we  require  a 
model  capable  of  accurately  predicting  spatially  distributed  hydraulic  variables  such  as  water 
depth  and  flow  velocity.  A  fundamental  constraint  to  the  development  of  such  a  model  is  the 
lack  of  suitable  calibration  and  validation  data  sets.  Whilst  it  is  now  possible  to  develop  high 
space/time  resolution  prediction  schemes,  current  data  capture  techniques  do  not  allow  the 
utility  of  such  developments  to  be  fully  assessed.  The  purpose  of  this  report  is  to  review 
recent  progress  in  hydraulic  modelling  and  examine  the  extent  to  which  recent  advances  in 
the  remote  sensing  of  flood  hydrology  have  the  potential  to  overcome  current  data  constraints 
thereby  leading  to  an  integration  of  remote  sensing  and  modelling  methodologies. 
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2.  REVIEW  OF  RECENT  PROGRESS  IN  HYDRAULIC  MODELLING 

2. 1  Background 

Progress  in  hydraulic  modelling  over  the  last  decade  has  led  to  considerable  improvements  in 
our  ability  to  simulate  river  flow  problems.  The  impetus  for  this  progress  has  come  from  a 
number  of  fields  and  incorporates  improvements  in  process  understanding,  mathematical  and 
numerical  developments  and  improvements  in  available  computational  power.  Large  scale 
flume  studies  (see  for  example  Knight  and  Shiono,  1993;  Sellin  and  Willetts,  1996)  have 
demonstrated  the  complexity  of  river  flow  processes  which  effectively  invalidate  one 
dimensional  models  for  certain  classes  of  flow  problem  (Knight  and  Shiono,  1996).  A  need 
has  therefore  arisen  to  develop  higher  resolution  hydraulic  models  (two  and  three 
dimensional)  to  account  for  the  impact  of  these  processes  during  river  flood  events. 
Concurrently,  developments  in  Computational  Fluid  Dynamics  (CFD)  have  led  to  the 
development  of  a  number  of  techniques  which  have  enabled  this  need  to  be  met.  Stable 
numerical  techniques  for  free  surface  flow  problems  have  been  developed.  For  example,  the 
Streamline  Upwind  Petrov  Galerkin  developed  by  Brookes  and  Hughes  (1982)  eliminates 
spurious  oscillations  from  numerical  solutions  of  the  Navier-Stokes  equations  without  adding 
artificial  diffiusivity.  Element-by-element  matrix  storage  techniques  (Carey  and  Jiang,  1986; 
Hervouet,  1992;  Binley  and  Beven,  1993)  have  reduced  computer  storage  requirements  and 
improved  efficiency,  allowing  spatially  complex  higher  dimensional  problems  to  be 
simulated.  Lastly,  and  of  critical  importance  to  flood  inundation  problems,  algorithms  have 
been  developed  to  simulate  dynamically  moving  inundation  extent  boundaries  (Lynch  and 
Gray,  1980;  King  and  Roig,  1988).  Thus,  the  last  decade  has  seen  a  significant  stimulus  to 
develop  two  dimensional,  or  higher,  numerical  models  for  flood  inundation  problems,  while 
at  the  same  time  a  series  of  technical  advances  has  rendered  this  goal  attainable. 

Two  and  three  dimensional  modelling  approaches  have  now  been  developed  for  application 
to  a  number  of  free  surface  flow  problems.  Such  models  solve  some  form  of  the  Navier 
Stokes  equations,  such  as  the  St.  Venant,  or  shallow  water,  equations  in  conjunction  with  a 
number  of  approximations  to  represent  such  processes  as  turbulence  and  boundary  fiiction 
(see  for  example  Hervouet  and  Van  Haren,  1996;  Lane,  1997).  An  appropriate  discretization 
of  the  domain  geometry  and  topography  is  developed  to  define  a  network  of  computational 
nodes  where  the  above  equations  can  be  solved  by  some  numerical  procedure.  Codes  have 
been  developed  using  a  variety  of  numerical  methods  ineluding  finite  difference  (Zeilke  and 
Urban,  1981),  finite  element  (Gee  et  al.,  1990;  Bates  et  al.,  1992)  and  finite  volume  (Lane  et 
al.,  1994).  This  solution  in  space  is  then  projected  forward  in  time  using  a  further  numerical 
procedure,  typically  a  finite  difference  scheme,  to  give  a  dynamic  simulation.  Given 
appropriate  boundary  conditions,  estimates  of  water  depth  and  flow  velocity  may  thus  be 
obtained  at  each  computational  node  at  each  time  step,  particular,  a  series  of  studies  have 
focused  on  the  development  of  hydraulic  models  for  long  reach  flood  inundation  problems. 
Until  recently,  application  of  hydraulic  models  to  free  surface  flow  problems  has  been 
characterised  by  three  specific  attributes:  the  scale  of  application  was  small  (0.5-2  km)  and 
focused  on  the  analysis  of  detailed  flow  patterns;  the  main  channel  was  rarely  resolved 
separately  from  the  floodplain  and  the  entire  domain  was  inundated  during  the  simulation. 
However,  most  floodplain  inundation  problems  of  practical  interest  occur  at  long  reach  scales 
(10-60  km)  and  involve  a  dynamic  flood  inundation  boundary.  As  a  consequence  there  is  a 
need  to  represent  dry  areas  within  the  computational  domain.  A  number  of  models  have 
therefore  been  developed  which  are  capable  of  application  at  this  scale  and  which  incorporate 
algorithms  to  account  for  wetting  and  drying  processes  (see  for  example  Gee  et  al.,  1990; 
Baird  et  al,  1992;  Feldhaus  et  al.,  1992;  LeClerc  et  al.,  1992;  Bates  et  al,  199^).  Given  the 
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scale  of  the  problem,  the  above  approaches  all  adopt  a  two  dimensional  finite  element 
approximation  to  the  flow  field  which  enables  complex  topography  to  be  represented  with  a 
minimum  number  of  computational  points. 

Such  simulations  are,  however,  by  no  means  straightforward  and  in  developing  an 
appropriate  flood  simulation  model  a  number  of  issues  have  to  be  considered.  Firstly,  the 
model  must  represent  a  complex  set  of  hydraulic  processes  operating  over  a  range  of  time  and 
length  scales.  These  range  from  turbulence  at  micro-scales,  for  example  1-10  cm  and  1-10 
seconds,  to  the  scale  of  the  flood  wave  itself,  which  may  be  hundreds  of  kilometres  in  length 
and  persist  for  several  days  (see  Bates  and  Anderson,  1993  for  a  more  complete  discussion). 
The  level  of  detail  at  which  processes  can  be  represented  is  dependent  on  the  scale  of  the 
simulation.  A  compromise  therefore  exists  between  processes  representation  and  scale,  with 
the  modeller  being  required  to  make  a  number  of  subjective  decisions  regarding  which 
processes  are  relevant  for  particular  applications.  Secondly,  field  observations  have  shown 
that  floodplains,  islands  and  channel  bars  possess  a  complex  micro  topography  that  may  have 
a  significant  effect  on  both  local  hydraulic  and  sediment  transport  processes  (Walling  et  al., 
1986;  Walling  et  al.,  1991)  as  well  the  inundation  extent  (Bates  and  Anderson,  1996)  for  a 
particular  recurrence  interval  flood.  Given  this  complexity  of  process  and  initial  condition 
specification,  hydraulic  models  must  simplify  actual  physical  processes.  This,  in  conjunction 
with  the  impossibility  of  obtaining  complete  a  priori  knowledge  of  parameters  such  as 
boundary  fiiction  in  distributed  models,  means  that  recourse  must  inevitably  be  made  to 
calibration  procedures  in  any  practical  application  to  achieve  a  correspondence  between 
model  predictions  and  field  data. 

2.2  Issues  of  Calibration  and  Validation 

Although  calibration  is  a  necessity,  it  may  lead  to  a  number  of  potentially  significant 
constraints  on  our  ability  to  further  develop  this  class  of  model.  In  effect  the  numerical 
procedure  attempts  to  solve  a  system  of  simultaneous  equations  where  there  are  many  more 
unknowns  than  equations.  The  objective  behind  estimating  these  unknowns  is  an  attempt  to 
obtain  a  match  with  some  test  data  set.  However,  given  that  data  availability  is  limited, 
particularly  when  compared  to  the  number  of  data  points  that  are  unknown,  the  degrees  of 
fi-eedom  present  within  the  modelling  procedure  are  such  that  the  solution  obtained  may  not 
be  unique.  Thus,  more  than  one  set  of  calibration  parameters  may  fit  the  available  data.  This 
is  particularly  true  of  inundation  extent  predictions.  Available  calibration  and  validation  data 
for  river  reaches  predominately  consist  of  observations  of  bulk  flow  properties  such  as 
discharge.  These  are  typically  available  only  at  the  reach  outflow,  as  in  the  River  Stour 
example  above.  It  may  be  entirely  possible  to  develop  a  number  of  calibrations  for  this 
problem  which  match  the  available  downstream  discharge  record  but  give  different 
predictions  of  flood  inundation  extent.  For  example.  Bates  et  al.  (1996)  studied  a 
hypothetical  channel/floodplain  domain  using  a  Monte  Carlo  procedure  to  generate  boundary 
fiiction  and  turbulent  viscosity  parameters  for  input  into  a  two  dimensional  finite  element 
model.  They  showed  that  a  range  of  calibrations  could  be  obtained  that  gave  equal 
acceptability  in  terms  of  their  ability  to  predict  the  downstream  outflow  hydrograph  but 
which  gave  inimdation  extents  varying  between  80  and  93%  of  the  total  floodplain  area. 
Thus,  incorrect  internal  process  representation  can  produce  a  'correct'  result  at  the  reach 
outflow  due  to  an  inability  to  constrain  the  calibration  procedure  to  a  unique  solution  (see  for 
example  Beven,  1989;  Konikow  and  Bredehoefl,  1992;  Fawcett  et  al.,  1995).  Moreover,  the 
results  of  this  study  show  that  we  are  unlikely  to  be  able  to  calibrate  complex  hydraulic 
models  using  hydrograph  data  alone.  Recent  modelling  developments  have  enabled  the 
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potential  for  prediction  to  far  outstrip  our  ability  to  collect  sufficient  calibration  and 
validation  data.  For  example,  the  River  Stour  simulation  discussed  above  generated 
approximately  1.17  x  lO’  items  of  data.  This  compares  to  an  available  data  set  of  248 
discharge  measurements.  Moreover,  these  data  are  not  spatially  distributed  and  we  are  thus 
unable  to  utilize  the  full  potential  of  the  model. 

A  lack  of  spatially  distributed  calibration  and  validation  data  therefore  currently  constrains 
our  ability  to  further  develop  high  resolution  flood  hydraulic  models.  In  order  to  determine  a 
research  design  for  future  model  development  we  require  sufficient  data  internal  to  the  model 
computational  domain  to  both  achieve  a  calibration  which  constrains  the  number  of  possible 
solutions  as  far  as  possible  and  to  validate  this  prediction  model  using  a  split-sample 
approach.  Due  to  the  problems  with  the  collection  of  hydraulic  data,  such  as  flow  velocity, 
during  flood  events  an  obvious  focus  for  attention  should  be  direct  observation  of  flood 
inundation  extent.  Not  only  is  this  variable  sensitive  to  small  changes  in  water  surface 
elevation  due  to  shallow  floodplain,  island  or  mud  flat  bed  gradients,  it  has  also  been  shown 
to  require  an  additional  set  of  approximations  over  and  above  those  needed  to  solve  the 
controlling  equations.  Its  accurate  prediction  is  therefore  a  good  test  of  model  capabilities 
and  of  significant  practical  interest.  In  this  context  remote  sensing  methods,  both  airborne 
and  satellite,  have  significant  potential  for  acquiring  such  information  at  a  spatial  scale 
commensurate  with  model  resolution. 

2.3  Key  Features  of  the  Study 

•  The  use  of  a  high  resolution  model,  in  both  space  and  time,  along  with  wetting  and 
drying  algorithms  for  representing  moving  flow  field  boundaries  allows  dynamic 
inundation  predictions. 

•  For  the  first  time  model  predictions  are  validated  on  this  scale  in  both  time  and  space 
using  multiple  data  sources.  The  data  sources  are  internal  stage  data  at  two  sites, 
supplied  by  the  USAGE  Missouri  River  District  (MRD),  and  satellite  imagery,  supplied 
by  the  Remote  Sensing  and  Geographic  Information  Systems  Center  at  USAGE  Cold 
Regions  Research  and  Engineering  Laboratory  (CRREL). 

•  The  influence  of  bathymetric  data  on  the  model  predictions  is  studied for  its  potential  in 
aiding  future  data  collection  strategy. 
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3.  A  TWO  DIMENSIONAL  MODEL  FOR  FLOODPLAIN  INUNDATION 

PREDICTION 

3. 1  2D  Model  Selection  for  the  Study 

The  TELEMAC  system  is  a  series  of  computer  programs  utilising  finite  element  techniques 
for  simulating  hydraulic  situations.  The  system  includes  pre-  and  post-processing  components 
and  offers  solutions  in  two-  and  three-dimensions.  It  also  contains  facilities  for  sediment  and 
contaminant  transport  (SUBIEF)  and  sand  transport  (TSEF).  The  full  system  is  further 
outlined  in  Figure  3-1.  All  the  components  of  the  system  have  common  file  formats  and  are 
written  in  the  high  level  language  FORTRAN  77.  This  enables  easy  file  transfer  between 
components  of  the  system  and  allows  model  users-to  modify  parts  of  the  code  as  desired  for 
specific  applications. 

The  TELEMAC  system  contains  both  two-  and  three-dimensional  versions.  The  three- 
dimensionaFversion  (TELEMAC-3D)  has  advantages  in  many  applications  where  vertical 
velocity  variations  are  important  such  as  small  scale  and  oceanic  studies  but  in  large  scale 
fluvial  applications  depth  averaged  calculations  are  adequate.  Indeed  the  equations  used  in 
nearly  all  fluvial  hydraulic  models  derive  fi'om  the  depth  averaged  equations  known  as  the 
Shallow  Water  or  St.  Venant  Equations,  TELEMAC-2D  being  no  exception.  The  third 
dimension,  though  initially  appealing,  would  simply  add  to  the  computational  and  data 
demands.  Descriptions  of  TELEMAC-3D  can  be  found  in  Hervouet  et  al.  (1994),  Hervouet 
and  Janin  (1994)  and  Hervouet  and  Van  Haren  (1996).  TELEMAC-2D  is  however  used 
throughout  the  work  presented  here  and  is  now  discussed  in  more  detail. 

Figure  3-1 :  The  TELEMAC-2D  modelling  system  (after  Hervouet  and  Lang,  1995). 


The  main  features  of  TELEMAC-2D  are  listed  by  Hervouet  et  al.  (1994)  as: 

•  structured  or  non  structured  meshing, 

•  use  of  Cartesian  or  spherical  co-ordinates, 

•  subcritical  and  supercritical  regimes  (with  hydraulic  jumps), 

•  various  momentum  source  terms  :  bottom  fiiction,  wind  stress  atmospheric  pressure, 
Coriolis  force, 

•  tiu-bulence  modelling  (k-epsilon  model), 

•  equation  on  temperature  or  a  substance  concentration, 

•  treatment  of  tidal  flats, 

•  many  types  of  boimdary  conditions  including  firee  slip  condition  and  incident  waves. 

Many  of  these  are  of  minimal  importance  for  fluvial  applications  but  are  necessary  for  the 
diverse  nature  of  cases  that  TELEMAC-2D  can  be  applied  to.  A  wide  range  of  test  cases  are 
illustrated  in  the  TELEMAC-2D  (version  3.0)  validation  document  (Cooper,  1996)  prepared 
to  substantiate  the  explicit  claims  made  about  the  applicability  and  accuracy  of  the  computer 
code.  The  test  cases  illustrated  range  from  the  Western  European  Coast  and  the  Mersey 
estuary  UK,  to  flows  around  bridge  piers  and  over  a  weir.  On  a  scale  more  relevant  to  this 
study  dam  breaks,  flow  at  a  river  confluence  and  the  simulation  of  a  flood  event  on  the  River 
Culm,  UK  are  also  shown.  A  new  version  of  TELEMAC-2D  is  released  annually.  All  of  the 
simulation  in  this  report  have  been  carried  out  using  version  3.0,  released  in  1995. 

3.2  Shallow  water  equations 
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Figure^-l  The  TELEMAC-2D  modelling  system  (after  Hervouet  and  Lang,  1995). 
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TELEMAC-2D  solves  the  shallow  water  equations  (SWE),  the  depth  averaged  version  of  the 
fully  three  dimensional  Navier  Stokes  equations  of  fluid  flow.  All  of  these  conditions  are 
commonly  met  in  rivers,  estuaries  and  seas  making  the  choice  of  these  equations  over  the  full 
three  dimensional  Navier  Stokes  equations  valid  for  the  applications  of  the  model  (Cooper, 
1996). 

The  choice  of  formulation  of  the  SWE  used  in  TELEMAC-2D  is  not  obvious.  A  conservative 
form  would  seem  better  but  divisions  by  the  water  depth  are  needed  to  produce  the  velocity 
field,  hence  causing  problems  in  dry  areas  of  the  model  domain.  Hence  a  non-conservative 
form  is  preferred.  Moreover  numerical  stability  analysis  also  favours  the  non-conservative 
version  of  the  equations.  Using  finite  element  methods  mass  conservation  can  be  ensured 
with  non-conservative  equations.  Two  versions  of  these  non-conservative  equations  have 
been  developed,  the  celerity-velocity  version  and  the  depth-velocity  version.  The  depth- 
velocity  version  are  marginally  favoured  as  they  are  easier  to  apply  mass-conservation 
techniques  to'(Hervouet  and  Janin,  1994).  These  equations  are  shown  below: 


—  +  u.grad(h)  +  h.div(u)  =  0 
(3 


Equation  3-1 


+  u.grad(u)  +  g-^  -  div(v.grad(u))  =  Sx  -  g-^ 
^  as.  as 


Equation  3-2 


dV  J/  \  ^  A-  !  C 

-—  +  u.grad(u)  +  g  —  -  div(v.grad(v))  =  Sy  -  g-— 
di  oy 


Equation  3-3 


where: 

h  :  water  depth 

u,  V  :  velocity  components 

g  :  gravity 

Zf  :  bottom  elevation. 

Sx,  Sy  :Source/sink  terms  (bottom  fnction,  wind,  etc.) 

V  :  eddy  viscosity 

h,  u  &  V  are  the  unknown  variables. 

Equation  3.1  being  the  continuity  equation,  3.2  and  3.3  the  force-momentum  equations.  Full 
derivations  can  be  foimd  in  Norton  et  al.  (1973). 

The  TELEMAC  system  uses  finite  element  methodology  to  solve  the  shallow  water  equations 
to  produce  values  of  water  depth  (h)  and  two  velocity  components  (u  and  v)  at  all  points  in 
the  model  domain.  To  achieve  this  the  domain  must  be  discretized  into  a  grid  or  mesh, 
usually  made  up  of  linear  triangles  (unstructured  grid)  for  flexibility.  The  mesh  is  created 
outside  TELEMAC  using  specialised  mesh  generation  software  such  as  SDRC  I-DEAS  or 
BALMAT.  The  mesh  is  then  integrated  with  the  topographic  data  to  produce  a  geometry  file 
that  is  then  used  in  the  TELEMAC-2D  simulation.  The  mesh  consists  of  a  series  of  node 
points,  the  vertices  of  the  triangles  with  a  elevation  (z)  value,  where  the  equatiojis  are  solved 
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and  the  triangles  themselves  are  the  elements.  The  geometry  notation  used  is  shown  in  Figure 
3-2. 


Figure  3-2  The  geometry  notation  used  in  TELEMAC-2D  (after  Hervouet  and  Van  Haren, 
1995).  _ 


3.2  Boundary  conditions 

There  are  two  main  types  of  boundary  condition  that  can  be  used,  solid  boundaries  and  liquid 
boundaries.  The  boundaries  described  here  are  round  the  sides  of  the  model  domain  and 
through  the  bed  of  the  model.  All  boundary  conditions  are  assigned  on  a  node  by  node  basis. 
Solid  boundaries  are  no  flux  (impermeable)  and  incorporate  a  friction  factor  (Hervouet  and 
Van  Haren,  1995).  Liquid  boundaries  allow  a  flux  across  them.  They  are  more  difficult  to 
deal  with  as  they  suppose  the  existence  of  a  fluid  area  that  is  not  part  of  the  calculation 
domain  but  can  however  influence  it.  This  influence  is  described  through  the  boundary 
condition.  There  are  four  types  of  liquid  boundary,  entry  and  exit  with  supercritical  flow 
(Froude  number  >1)  and  entry  and  exit  with  subcritical  flow  (Froude  number<l).  Incident 
waves  and  prescribed  flowrates  can  be  incorporated  through  these  4  boundary  types. 
Hervouet  and  Van  Haren  (1995)  describe  these  boundary  conditions  more  fully. 

3.4  Physical  parameter  options 

Physical  representation  of  parameters  is  essential  in  models  such  as  TELEMAC-2D  in  order 
to  apply  them  to  different  applications.  TELEMAC-2D  includes  numerous  physical 
parameters.  The  most  important  are  discussed  in  this  section. 

There  are  six  options  in  TELEMAC-2D  for  representing  bed  friction,  these  are: 

•  No  friction 

•  Linear  friction 

•  Chezy's  law 

•  Strickler's  law 

•  Manning's  law 

•  Nikuradse's  law 
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The  applied  friction  coefficient  is  converted  to  force  terms  in  the  x  and  y  directions  at  each 
computational  node  which  is  then  fed  into  the  Shallow  Water  Equations  (in  the  source  terms 
Sx  and  Sy  -  see  section  3.2,  Equation  3-2  and  Equation  3-3).  In  reality  the  bed  friction  force 
is  a  quadratic  function  of  velocity  so  the  no  friction  and  linear  law  are  rarely  used  in  practical 
applications  (Hervouet  and  Van  Haren,  1996).  The  Chezy,  Strickler  and  Manning  laws  are  all 
closely  related  and  utilise  the  quadratic  function  mentioned  and  are  all  described  more  fully 
by  Hervouet  and  Van  Haren  (1995).  Nikuradse's  law  calculates  a  Chezy  coefficient  from  the 
water  depth  (h)  and  grain  size  of  the  bed  material  which  is  then  converted  into  a  force  term. 
Which  of  these  laws  that  is  used  is  not  very  important  as  they  are  very  closely  related  and 
friction  coefficients  can  easily  be  converted  between  them.  The  choice  rests  with  the  model 
user. 

Applications  of  TELEMAC-2D  to  river  and  estuaries  generally  involve  wetting  and  drying 
areas  of  the  model  domain.  For  example  the  tidal  simulation  in  the  Mersey  estuary  and  the 
flood  event  *on  the  River  Culm,  both  described  by  Cooper  (1996),  involve  this  type  of 
behaviour.  The  ability  of  the  model  to  deal  with  this  sort  of  behaviour  is  therefore  of  vital 
importance  but  is  problematic.  The  behaviour  in  the  dry  zones,  where  divisions  by  the  water 
depth  (h)  in  the  calculations,  can  cause  spurious  terms  to  appear  as  h  tends  towards  zero. 

Two  solutions  to  this  problem  are  available  in  TELEMAC-2D,  both  described  in  more  detail 
by  Hervouet  and  Van  Haren  (1995), 

•  solving  the  equations  ever3where  and  coping  with  the  spurious  terms. 

•  removing  the  dry  zones  from  the  computational  domain. 

The  first  is  the  simplest  but  corrections  must  be  applied  in  the  dry  zone.  In  dry  areas  the  water 
surface  gradient  becomes  that  of  the  bottom  topography  but  this  cannot  be  allowed  to  act  as  a 
driving  force  in  the  momentum  equation. 

The  second  method  removes  the  dry  zone  from  the  computational  domain  and  is  often  called 
the  "moving  boundary  technique".  In  TELEMAC-2D  this  is  achieved  efficiently,  avoiding 
the  need  to  redefine  the  mesh  at  each  time  step,  by  keeping  the  elements  in  the  mesh  but 
cancelling  their  existence  through  the  use  of  an  array  set  to  0  for  dry  elements  and  1  for  all 
others. 

Partially  dry  elements  are  another  important  area,  especially  where  the  elements  are  large 
such  as  in  the  Missouri  River  model.  These  are  coped  with  in  TELEMAC-2D  by  a 
sophisticated  method  that  allows  the  water  depth  to  go  to  zero  within  an  element.  This 
compares  favourably  with  the  usual  techniques  of  keeping  the  elements  fully  wet  or 
excluding  partially  dry  elements  (Figure  3-3). 

3.5  Summary  of  the  TELEMAC-2D  system 

TELEMAC-2D  is  a  high  resolution  space/time  distributed  hydraulic  model  that  solves  the 
Shallow  Water  Equations  for  fluid  flow  using  finite  methodology.  The  model  can  be  used  in  a 
wide  variety  of  scenarios  including  those  involving  wetting  and  drying  fronts  within  the 
model  domain.  Bed  friction  and  turbulence  are  represented  in  the  model  through  the  use  of 
physically  based  parameters. 

TELEMAC-2D  can  therefore  be  seen  to  be  well  specified  for  the  type  of  fluvial  application 
that  is  of  interest  in  this  study.  The  code  has  been  shown  to  work,  through*  the  report  of 
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Cooper  (1996),  in  a  wide  range  of  cases.  Success  in  any  individual  situation  is  however 
dependant  on  the  data  provided,  physical  and  numerical  parameters  used.  How  these 
important  factors  have  been  determined  in  the  Missouri  River  model  case  is  the  subject  of  the 
next  section.  _ 


(a) 


Free  surface 
elevation 


Figure  3-3  Partially  dry  elements:  (a)  the  real  situation,  (b)  element  remaining  fully  wet 
introducing  spurious  momentum  terms,  (c)  element  exclusion  method  and  (d)  TELEMAC-2D 
method  keeping  the  element  partially  dry  (after  Price,  1997). _ _ 


A 
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4.  THE  MISSOURI  RIVER  MODEL 

4. 1  Missouri  River  Reach 

The  reach  of  the  Missouri  River  being  used  for  this  modelling  study  is  that  from  Gavins  Point 
Dam,  South  Dakota,  to  the  gauging  station  at  Maskell,  Nebraska  (Figure  4-1).  The  reach 
covers  river  miles  81 1  to  776  making  it  35  miles  or  55  km  long.  The  channel  varies  in  width 
between  300m  and  1200m.  The  channel  slope  is  very  low,  dropping  only  about  12  metres 
along  its  55  km  length,  giving  a  gradient  of  0.02%.  The  bed  material  in  this  channel  is  sand 
which  is  fairly  mobile  but  the  channel  banks  have  been  strengthened  or  stabilised  along  much 
of  the  reach.  There  are  several  islands  along  the  reach  several  of  which  are  permanent.  There 
is  one  major  tributary,  the  James  River,  that  joins  the  main  stem  at  river  mile  800  adjacent  to 
the  James  River  Island. 

The  flow  out  of  Gavins  Point  Dam  is  regulated  to  minimise  the  risk  of  flooding  downstream, 
hence  the  reach  being  modelled  is  very  unlikely  to  attain  out  of  bank  conditions.  The  flow  in 
the  James  River  is  however  naturally  variable.  The  model  of  this  system  can  therefore  remain 
as  channel  only.  River  flow  data  is  available  at  several  points  along  the  reach  on  an  hourly 
basis.  Table  4-1  shows  the  data  availability  and  gauge  locations  (also  Figure  4-1). 


Table  4-1  Gauge  station  location  and  data  availability. 


Gauging  Station 

Data  Available 

Gavins  Point  Dam,  SD. 

Yankton,  SD. 

Gayville,  SD. 

Maskell,  Neb. 

Scotland,  SD. 
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Figure  4-la)  Location  and  (b)  detail  of  the  study  reach  of  the  Missouri  River,  USA.  The 
modelled  reach  is  from  Gavins  Point  Dam  to  Maskell  and  includes  the  gauging  stations  at 
Yankton  and  Gayville. 
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The  flow  from  the  James  River  gauge  at  Scotland,  South  Dakota,  is  routed  using  a  simple 
one-dimensional  kinematic  wave  model  to  the  required  location  on  the  main  model  input 
boundary  when  necessary. 

The  observed  data  for  the  Missouri  model  is  of  very  high  quality  and  quantity.  The  available 
data  sets  are  of  special  value  because  they  allow  both  detailed  point  comparisons  using  the 
stage  data  at  two  sites,  Yankton  and  Gayville,  and  spatially  distributed  comparisons  using  the 
satellite  imagery.  Using  this  data  the  performance  of  the  Missouri  River  model  can  be 
comprehensively  analysed. 

4.2  Producing  the  finite  element  mesh 

The  finite  element  mesh  for  the  model  must  be  created  by  first  defining  the  boundary  of  the 
model  domain  and  secondly  discretizing  this  area  into  elements. 

The  boundary  of  the  2D  model  was  defined  by  digitising  round  the  edge  of  the  river  as 
represented  on  United  States  Department  of  the  Interior  7.5  minute  quadrangle  series  maps. 
The  use  of  wetting  and  drying  algorithms  in  the  model  enable  the  flow  field  boundary  to  be 
calculated  within  this  outer  boundary  thus  allowing  this  boundary  to  be  fairly  loosely  defined. 
The  digitised  boundary  was  then  converted  to  Universal  Transverse  Mercator  (UTM)  metre 
co-ordinates  for  compatibility  with  the  metric  scale  used  by  TELEMAC-2D.  The  boundary  of 
the  model  (with  TELEMAC-2D’s  metric  scale),  which  includes  three  permanent  islands,  is 
shown  in  Figure  4-2. 


A  mesh  of  triangular  elements  must  be  made  within  the  model  domain  to  facilitate  the  finite 
element  solution  technique.  The  finite  element  mesh  used  with  the  Missouri  model  were 
generated  inside  the  prescribed  boundary  using  the  mesh  generation  package  BALMAT  and 
consisted  of  near  equilateral  triangles  in  order  to  increase  the  accuracy  and  minimise  mass 
conservation  errors.  This  consisted  of  5969  nodes  and  10437  elements  with  an  average  size 
of  approximately  100  m^  and  is  shown  in  Figure  4-3. 

Figure  4-3  Finite  element  mesh  of  the  Missouri  River  from  Gavins  Point  Dam  to  Maskell. 
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4.3  Bathymetric  Data 

The  bathymetric  data  in  a  2D  hydraulic  model  is  one  of  the  most  important  factors  in 
attaining  high  quality  simulations.  The  quantity  and  quality  of  the  bathymetric  data  is 
therefore  of  utmost  significance. 

Along  this  reach  of  the  Missouri  River  the  bathymetry  was  obtained  by  the  MRD  in  1995  by 
echo  sounder  surveys  of  channel  cross  sections.  In  total  343  cross  sections  were  available 
along  the  55  km  reach  creating  an  average  spacing  of  just  over  150m.  Each  cross  section  was 
made  up  of  between  30  and  200  data  points,  depending  on  section  length.  The  typical  spacing 
between  data  points  along  a  section  was  6m.  All  elevation  data  was  converted  to  metres. 
Bathymetric  data  was  missing  on  the  cross  sections  where  the  bed  was  above  the  water 
surface,  such  as  on  sand  banks  and  mud  flats.  In  these  regions  the  topographic  elevation  was 
estimated  from  United  States  Department  of  the  Interior  7.5  minute  quadrangle  maps,  except 
on  the  permanent  island  in  the  model  where  no  such  data  is  required. 

The  bathymetric  data  must  now  be  applied  to  the  finite  element  mesh.  This  is  the  process  of 
interpolating  between  topographic  data  points  and  assigning  elevation  (z)  values  to  the  nodal 
points  in  the  mesh.  This  produces  a  geometry  file  for  the  model  simulation.  The  interpolation 
is  usually  carried  out  using  STBTEL,  the  TELEMAC  sub-program.  This  uses  a  quad- 
directional  interpolation  routine,  taking  the  nearest  data  point  in  the  four  quadrants  around 
each  node,  weighting  for  distance  from  the  node  and  combining  to  give  the  nodal  elevation. 
An  alternative  interpolation  routine  has  been  written  at  Bristol  for  use  with  cross  sectional 
data  of  river  channels.  This  is  a  linear  interpolation  down  the  line  of  greatest  depth  between 
sections  and  improves  the  definition  of  this  line  in  such  cases  as  the  Missouri  River.  A 
portion  of  the  models  topography  is  shown  in  Figure  4-4. 
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Figure  4-4  An  example  portion  of  the  models  bathymetric  representation.  This  section  is 
immeadiately  upstream  of  Yankton. _ _ 
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4.4  Boundary  Condition  Specification 

The  boundary  conditions  of  the  model  are  very  important  to  the  simulations  run.  There  are 
several  boundaries  to  each  model  that  must  be  set  with  care  prior  to  the  model  run.  These 
boundaries  are  divided  into  two  types,  liquid  boundaries  and  solid  boundaries. 

With  the  Missouri  model  the  liquid  boundaries,  allowing  flow  across  them,  are  at  the  top  and 
bottom  ends  of  the  reach  and  on  the  James  River  tributary  inflow.  Model  inflows  are 
prescribed  as  flowrates,  at  Gavins  Point  Dam  and  from  the  James  River,  and  outflows  as 
water  surface  elevations,  at  Maskell.  This  produces  a  well  posed  problem  for  fluvial  flows 
according  to  the  theory  of  characteristics  (Hervouet  and  Van  Haren,  1995).  All  measured 
flowrates  were  converted  to  m^/s  (cumecs)  and  stage  values  to  metres. 

Solid  boundaries  allow  no  flux  across  them.  They  are  found  down  the  sides  of  the  reach  and 
through  the  bed.  The  side  boundaries  in  the  Missouri  model  are  set  as  slip  boundaries, 
allowing  a  velocity  along  them.  This  is  justified  by  the  size  of  the  elements  making  flow 
overly  restrained  in  areas  of  channel  constriction  when  the  more  realistic  no-slip  boundaries 
are  imposed.  The  flexible  boundary  of  the  flow  field  removes  this  problem  for  large  portions 
of  the  reach  where  such  an  argument  is  irrelevant. 

4.5  Physical  Parameter  Specification 

The  two  most  important  physical  parameters  in  fluvial  hydraulic  models  are  generally  agreed 
to  be  bed  fiiction  and  turbulence  (Baird  and  Anderson,  1990;  Bates  et  al,  1992).  The  theory 
of  the  two  in  TELEMAC-2D  has  been  discussed  in  section  3.4.  They  are  the  only  two 
considered  in  any  detail  here. 

Given  the  channel  only  nature  of  the  model  and  lack  of  additional  information  the  bed 
fiiction  and  turbulence  parameters  were  defined  as  constant  throughout  the  entire  reach. 
These  are  obviously  simplifications  but  should  be  adequate  as  a  first  approximation.  The  bed 
fiiction  is  always  defined  using  Manning's  law  enabling  the  standard  Manning's  'n'  measure 
of  flow  resistance  to  be  used.  The  turbulence  was  defined  using  the  zero-equation  velocity 
diffusivity  representation. 
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5.  COMPARISON  OF  MODEL  PREDICTIONS  TO  OBSERVED  DATA 

5. 1  Methods  of  Comparison 

The  data  available  for  this  study  is  stage  data  at  two  sites  within  the  model  domain  which  can 
be  considered  internal  but  1-D  and  numerical.  The  satellite  images  are  also  internal  but  2-D 
and  numerical  (except  for  example  of  process  validation  that  shall  be  shown  later).  The 
satellite  imagery  is  therefore  the  more  powerful  validation  data  type.  How  well  the  observed 
and  predicted  match  shall  be  shown  in  the  following  sections  and  at  the  end  this  idea  of 
strength  of  validation  shall  be  revisited  in  the  light  of  the  results.  Firstly  however  the  methods 
of  model  application  used  in  this  study  and  elsewhere  are  reviewed. 

Usually  when  hydrological  models  are  applied  to  actual  scenarios  a  two  stage  process  of 
model  calibration  and  validation  is  carried  out.  Calibration  is  the  process  by  which  parameter 
values  are  varied  within  reasonable  ranges  until  the  differences  between  observed  and 
computed  values  are  minimised  (Konikow  and  Bredehoeft,  1992).  Theoretically  speaking, 
physically  based  model  should  not  need  calibrating.  The  estimation  and  application  of 
physically  realistic  parameter  values  for  the  model  should  enable  the  model  to  perform  well 
without  any  further  manipulation.  However  the  numerous  simplifications  involved  in 
modelling  and  problems  in  parameter  estimation  mean  that  such  a  notion  is  unrealistic. 
Calibration  is  therefore  performed  on  virtually  all  physically  based  hydrologic  models. 
Following  the  calibration  phase  the  model  must  then  be  validated.  In  most  cases  the 
calibrated  parameter  values  are  used  again  on  a  different  portion  of  the  flow  record  to  see  if 
they  enable  a  good  match  between  observed  and  predicted  variables.  If  they  do  then  the 
model  and  parameters  are  considered  validated  and  can  be  used  with  confidence  for 
prediction  of  future  events.  If  not  then  calibration  and  possibly  some  stage  of  model 
formulation  must  be  repeated  until  the  model  is  validated.  It  should  be  noted  at  this  point  that 
validation  is  used  here  to  mean  that  a  reasonable  comparison  is  obtained  between  observed 
and  predicted  data  rather  than  that  the  model  has  been  shown  to  be  a  accurate  representation 
of  the  system. 

With  this  application  of  TELEMAC-2D  the  highly  structured  methodology  for  calibration 
and  validation  outlined  above  is  not  entirely  appropriate  for  several  reasons: 

•  this  study  is  for  research  of  model  behaviour, 

•  the  niunber  and  variety  of  data  sources  make  its  application  difficult, 

•  there  is  no  need  to  produce  a  model  capable  of  predicting  future  events  at  this  stage. 
The  two  processes  are  therefore  not  used  as  distinct  entities  in  the  following  investigation. 
Instead  a  hybrid  of  the  two  forms  was  used  enabling  all  benefits  of  the  calibration-validation 
procedure  to  be  accrued  but  also  much  more.  The  model  was  therefore  run  using  a  wide  range 
of  bed  friction  values  over  several  sections  of  the  flow  record.  This  enabled  the  calibrated 
parameter  values  to  be  found  for  each  comparison  data  set  (2  internal  stage  gauges  and 
possibly  some  areas  on  a  satellite  image).  Cross  comparison  of  parameter  values  and  the 
assessment  of  errors  when  the  model  was  calibrated  onto  one  data  source  allowed  an 
assessment  of  the  model/parameter  performance  more  complete  and  powerful  than  with  the 
traditional  calibration-validation  procedure. 

The  time  period  being  simulated  for  this  comparison  is  around  6th  June  1994  synchronous  to 
the  LANDSAT-TM  satellite  image  of  the  area  that  will  be  used  in  section  5.3.  The  flow 
record  around  this  date  showed  very  little  (<2%)  variation  in  the  data  set.  The  simulation  was 
therefore  taken  as  a  simple  steady  state  computation,  simplifying  the  analysis  of  results.  All 
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the  following  comparisons  are  therefore  done  on  a  single  point  in  time  basis,  e.g.  a  single 
water  level  for  each  point  in  the  domain  for  each  parameter  set. 

5.2  Comparison  of  model  prediction  with  internal  stage  data 

Two  internal  gauges  for  water  level  are  available  for  the  reach  of  the  Missouri  River 
described  in  section  1.  The  gauges  are  at  Yankton,  7.5  km  down  the  reach,  and  Gayville,  21.5 
km  down  the  reach  (Figure  4-1).  Hourly  stage  records  are  available  for  these  gauges. 

The  sensitivity  of  the  model  prediction  of  the  water  depth  at  these  two  locations  has  been 
shown  to  be  dependant  predominantly  on  the  bed  friction  parameter.  This  is  the  only 
parameter  varied  in  the  following  analysis.  Other  parameters  would  have  only  minimal  effect 
but  increase  complexity  greatly.  The  sensitivity  of  water  depth  to  friction  is  virtually  linear  at 
both  interna]  gauge  sites  (Figure  5-1)  and  in  the  expected  direction,  i.e.  water  depth 
increasing  as  friction  increases. 


This  finding  allows  the  use  of  McCuen's  (1973)  linearized  sensitivity  equation  but  a  simple 
regression  equation  is  more  useful  for  applying  the  results  in  this  analysis  and  is  therefore 
used  here.  Equations  were  produced  for  the  model  behaviour  at  both  sites,  Yankton  and 
Gayville,  relating  the  error  in  water  level  prediction  to  the  bed  fnction  (Manning's  'n') 
parameter.  The  equations,  calculated  using  the  statistical  computing  package  MINITAB,  can 
be  written  thus: 

Yankton  Error  =  -0.410  +  34.72  n 
Gayville  Error  =  -0.976  +  31.93  n 

Alternatively  they  can  relate  the  water  depth  to  the  bed  fnction  thus: 

Yankton  Water  Depth  =  4.87  +  34.72  n 
Gayville  Water  Depth  =  2.91  +  31.93  n 

The  confidence  in  these  regression  tines  is  huge  with  tiny  residuals  (<0.01 5m)  in  all  cases. 
These  equations  can  therefore  be  used  as  a  very  powerful  tool  for  estifhating  model 
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predictions  for  'n'  values  at  which  the  model  has  not  been  run,  although  it  should  not  be 
extrapolated  beyond  the  range  of  'n'  values  used  to  make  the  relationship  (i.e.  0.01  to  0.04). 
At  present  this  is  only  shown  at  two  points  but  could  perhaps  be  applied  at  every  node. 

It  can  easily  be  calculated  from  the  above  equations  that  the  'n'  value  to  eliminate  the  error  in 
water  surface  elevation  at  Yankton  is  0.0118  whereas  at  Gayville  it  is  0.0306.  The  power  of 
the  above  regression  equations  can  now  be  utilised  to  calculate  the  error  in  the  other 
observation  whilst  one  is  correct.  Whilst  the  stage  at  Yankton  is  predicted  correctly  at 
Gayville  the  prediction  is  0.60m  too  low.  Reversing  this,  whilst  Gayville  is  predicted 
correctly  the  prediction  at  Yankton  is  0.65m  too  high.  In  percentage  terms  these  errors  are 
18.26%  and  10.96%  of  the  expected  water  depth  at  the  two  points  respectively. 

Taking  the  value  central  to  the  above  two  estimates  of  Manning's  'n'  should  enable  the 
approximate  calculation  of  the  joint  minimum  errors  in  the  predictions.  This  value  of  'n'  is 
0.0212  and  'produces  an  error  at  Yankton  of  0.33m  and  at  Gayville  of  -0.30m.  Slight 
adjustment  of  'n'  could  equalise  these  errors  but  would  be  irrelevant.  These  errors  can  be 
expressed  again  as  percentages  of  the  water  depth  such  that  the  error  at  Yankton  is  5.89%  and 
at  Gayville  is  8.37%.  Figure  5.2  shows  a  downstream  section  of  the  comparison  between 
observed  and  predicted  water  surface  values. 


Interestingly  the  value  of  rate  of  change  of  the  water  surface  (error  term  or  water  depth)  is 
different  at  the  two  gauging  stations.  Thus  a  fixed  change  in  the  bed  fnction  produces  a 
different  change  in  the  water  surface  elevation  at  the  two  stations.  Investigation  into  how  this 
type  of  behaviour  occurs  over  the  whole  reach  would  be  useful  in  determining  more  about 
model  behaviour  or  what  causes  this  phenomenon. 

The  comparison  between  observed  an  predicted  stage  values  has  shown  that  the  model  is 
capable  of  predicting  this  variable  to  about  0.3m  at  two  sites  simultaneously.  This  is  a  good 
result  given  the  simplicity  of  the  model  set  up.  Given  further  work  on  the  parameterization 
this  could  no  doubt  be  reduced  significantly. 

5.3  Comparison  of  model  inundation  predictions  with  satellite  imagery 
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Many  types  of  satellite  imagery  have  been  used  widely  in  studies  of  fluvial  and  coastal 
hydraulics.  The  wide  coverage,  ease  of  use  and  wealth  of  information  contained  in  the  images 
makes  them  ideal  for  many  uses  in  this  field.  Redfem  and  Williams  (1996)  review  the 
available  sensors  and  their  applications  in  (primarily  coastal)  hydraulics.  Previously  a  lot  of 
work  has  been  carried  out  on  the  remote  sensing  of  flood  extent  because  of  the  immense 
hazard  and  mitigation  costs  on  floodplain  developments.  Rango  and  Salomonson  (1974) 
show  that  the  areal  extent  of  flooding  can  be  mapped  using  near-infrared  sensors  on  ERTS  1 
(Earth  Resources  Technology  Satellite  -  later  renamed  Landsat).  More  recently  Imhoff  et  al. 
(1987)  use  SAR  (Synthetic  Aperture  Radar)  and  Landsat  MSS  (Multispectral  Scarmer)  for 
mapping  flood  extent  and  damage  in  Bangladesh.  SAR  has  been  shown  to  very  accurate  for 
flood  boundary  delineation  by  Biggin  and  Blythe  (1996)  on  the  River  Thames  in  the  UK. 
Satellite  remote  sensing  of  the  1993  floods  on  the  Mississippi  and  Missouri  Rivers  in  the 
USA  with  both  SAR  (Brackenridge  et  al,  1994)  and  Landsat  TM  (Thematic  Mapper)  have 
further  shown  the  potential  of  such  methodology. 

Satellite  data  has  also  been  used  to  measure  flood  stages  as  illustrated  by  Koblinsky  et  al 
(1993)  using  the  U.S.  Navy’s  Geosat  radar  altimeter  on  the  Amazon  basin  and  Brackenridge 
et  al  (1994)  using  SAR  data  and  topographic  maps  on  the  Mississippi.  Both  methods  have 
quite  large  errors  associated  with  them  at  present. 

Plumes  of  sediment  rich  water  can  also  be  identified  as  they  have  a  higher  reflectance  than 
clear  water  in  the  visible  region  of  the  spectrum  (Lillesand  and  Kiefer  1987).  Brackenridge  et 
al  (1994)  use  this  factor  to  highlight  levee  breaches  along  the  Mississippi  River  during  the 
1993  floods,  the  breaches  acting  as  a  sediment  sources,  and  explain  sediment  deposition 
patterns.  Other  forms  of  pollution  and  thermal  emissions  can  also  be  traced  (Lillesand  and 
Kiefer,  1987). 

The  availability  of  satellite  imagery  for  the  modelled  reaeh  of  the  Missouri  River 
synchronous  to  the  flow  records  has  enabled  several  fundamental  research  objectives  to  be 
addressed.  Data  from  satellite  images  is  ideal  for  use  with  2-D  fluvial  hydraulic  models  being 
of  a  scale  commensurate  to  the  model  resolution  and  being  widely  spatially  distributed.  In 
this  study  such  data  has  the  potential  to  be  used  for: 

•  validation  of  inundation  extent  predictions, 

•  identify  regions  of  error, 

•  calibration  data, 

•  process  validation. 

Remotely  sensed  data  has  the  potential  to  be  used  more  widely  in  hydraulic  modelling  for 
topographic  and  physical  parameterization  (Bates  et  al,  1997)  but  these  are  ongoing  research 
themes. 

Two  sites  at  which  to  analyse  model  immdation  predictions  were  chosen.  These  were  selected 
specifically  because  of  their  complex  bathymetry  and  topography,  involving  mud  flats, 
permanent  and  temporary  islands.  The  location  of  the  two  sites  are  shown  on  Figure  5-3.  Side 
by  side  comparisons  of  the  observed  and  predicted  inundation  are  to  be  made  using  different 
fHction  values  at  both  sites  followed  by  an  overlay  showing  the  geo-referencing  that  can  be 
done  to  the  results. 
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Figure  5-3:  Locations  of  the  two  sites  to  be  used  in  the  comparison  of  model  predictions  with 
satellite  imagery. 


The  side  by  side  comparisons  of  the  model  predictions  against  the  Landsat  TM  image,  using 
bed  friction  values  of  0.01,  0.02,  0.03  and  0.04,  are  shown  in  Figure  5-4  and  Figure  5-5  for 
the  two  sites.  These  figures  plot  the  satellite  image  against  the  predicted  flow  field  boundary 
as  this  is  the  model  result  that  can  be  directly  compared  to  the  image.  From  these  it  can  be 
seen  that  the  model  predictions  are  generally  good,  producing  a  close  spatial  match  against 
the  observed  data.  This  allows  calibration  of  the  model  to  be  carried  out  with  such  data  by 
varying  bed  friction  until  there  is  a  close  spatial  match  between  observed  and  predicted.  This 
variable  behaviour  is  caused  by  the  gradient  of  the  topography  around  the  water  surface.  If 
the  topography  has  a  shallow  gradient  then  changes  in  the  bed  friction,  changing  the  water 
depth,  causes  changes  in  the  inundated  areas.  Where  the  topographic  gradient  is  steeper  then 
the  same  changes  in  water  depth  shall  have  virtually  no  impact  on  inundated  area.  A  problem 
arises  in  this  specific  case,  as  far  as  calibrating  the  model  on  this  data,  because  the 
topography  above  the  water  surface  was  only  estimated  from  US  Department  of  the  Interior 
maps  and  is  therefore  probably  not  of  the  accuracy  required  to  enable  much  confidence  to  be 
placed  in  the  predictions  of  immdated  area  on  islands. 

Figure  5-4:  Comparison  of  Landsat  TM  image  and  4  alternative  model  predictions  with 
different  bed  friction  (Manning’s  ‘n’)  values  at  site  1. 

Figure  5-5:  Comparison  of  Landsat  TM  image  and  4  alternative  model  predictions  with 
different  bed  friction  (Manning’s  ‘n’)  values  at  site  2. 

By  geo-referencing  the  model  results  and  satellite  image  an  overlay  of  the  two  can  be  made. 
An  example  of  such  is  shown  in  Figure  5-6  for  a  section  of  this  reach.  This  detailed 
comparison  between  observed  and  predicted  flow  field  boundary  enables  problem  areas  to  be 
pinpointed  and  further  model  development  to  proceed  by  further  data  collection  or 
adjustments  to  the  model  in  these  areas. 
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Figure  5-4*  Comparison  of  Landsat  TM  image  and  4  alternaUve  model  predictions  with 
different  bed  friction  (Manning’s  ‘n’)  values  at  site  1  for  event  1. 


I'lgurc  5-6  An  overlay  of  the  model  prediction  of  flow  field  boundary  onto  the  Landsat  TM 
l^nnayc  lor  a  rjpgion  around  Goat  Island  (including  sites  5  and  6)  for  e\  cnt  1. 

.‘'.-r  SuDlDUirV 

1  hesc  Hrsl  simulated  (low  conditions  ha\^e  shown  the  both  the  excellent  predictions  of  the  2- 
D  h>'draulic  model  and  the  utility  of  the  data  sources  to  validate  its  performance.  The  internal 
^lage  data  was  matched  at  single  points  but  not  so  closely  at  both  points  simultaneously.  The 
inundation  extent  was  validated  against  the  satellite  images  very  well.  Some  areas  showed  up 
pioblems  in  the  models  bathymetry/topography.  Given  possible  future  refinements  in  the 
model  set-up  the  internal  stage  data  should  be  able  to  be  matched  simultaneously  at  both  sites 
piiite  easily.  Matching  all  the  How  field  boundaries  shall  be  more  difficult  given  tlieir 
oependence  on  the  model  bathymetry.  This  is  however  expected  giN  Cii  tlie  relaii\'e  strength  of 
the  two  validation  tvpes  as  discussed  earlier. 
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6.  DERIVING  BATHYMETRIC  DATA  REQUIREMENTS 

6. 1  Initial  considerations 

The  high  cost  of  collecting  bathymetric  and  topographic  data  must  be  justified  by  attaining 
high  quality  model  results  that  could  not  be  accomplished  with  a  more  basic  representation. 
The  impact  of  the  density  of  cross  section  on  the  models  bathymetry  and  model  predictions 
are  investigated  in  this  section. 

Four  realisations  of  the  model  bathymetry  were  created  using: 

•  all  the  available  cross  sections  (150  m  spacing), 

•  alternate  cross  sections  (300  m  spacing), 

•  one  in  four  cross  sections  (600  m  spacing), 

•  one  in  eight  cross  sections  (1200  m  spacing). 

From  each  of  these  data  sets  the  model  bathymetry  was  created  at  the  node  points  on  the  low 
resolution  nfesh  using  the  directional  interpolator  described  in  section  4.3.  The  model  was 
then  run  once  with  each  of  these  representations  for  the  conditions  of  event  1  with  a  bed 
friction  value  of  0.025. 

The  impact  of  these  changes  on  the  model  bathymetry  was  profound.  Figure  6-1  shows  the 
impacts  in  one  example  area  (site  6).  As  can  be  clearly  seen  the  detail  of  the  bathymetric 
representation  in  the  model  degrades  rapidly  as  the  number  of  cross  sections  decrease. 
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Figure  6-1  Impact  of  degrading  the  number  of  cross  sections  used  on  the  model  bathymetry. 
The  spacing  degrades  from  150  m  to  1200  m  across  the  range  illustrated. _ _ 
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6.2  Comparison  to  observed  data 

The  influence  of  changing  the  level  of  bathymetric  data  on  the  water  surface  elevation 
predictions  at  Yankton  and  Gayville  is  shown  in  Table  6-1.  From  this  table  it  can  be  seen  that 
at  Yankton  the  level  of  data  provision  has  minimal  effects  of  the  predictions,  varying  them  by 
only  O.lm  compared  to  an  actual  error  of  approaching  0.5m.  At  Gayville,  however, 
decreasing  the  bathymetric  produces  a  continuous  increase  in  water  surface  elevation.  The 
magnitude  of  the  change  over  the  runs  is  0.72m,  spanning  the  observed  value.  The  half  data 
set  run  produces  a  prediction  of  water  surface  elevation  at  Gayville  within  0.015m  of  the 
observed  value. 


Table  6-1  Water  surface  elevations  at  Yankton  and  Gayville  for  simulations  with  the  various 


Location 

Full 

Half 

Quarter 

Eighth 

Observed 

Yankton 

■5352.667  5'. 

V352l721-j;::.:, 

352:i96  i?. 

Gayville 

B49;854  . 

:349.990>::i,-i> 

B503;98'f®^ 

35Q.032 

From  these  results  there  appears  to  be  little  benefit  in  using  the  highest  density  data  set  as  the 
results  are  equally  good  using  half  the  data.  The  quarter  and  eighth  runs  do  however  show 
decay  in  the  quality  of  the  results  at  Gayville.  The  more  important  results  are  however  the 
spatially  distributed  ones  that  are  now  investigated. 

Once  more  the  side  by  side  comparisons  have  been  made  for  model  predictions,  this  time 
with  the  4  levels  of  bathymetric  data,  against  the  Landsat  TM  image  at  the  two  sites  of 
interest  (Figure  6-2  and  Figure  6-3).  In  each  case  the  full  bathymetry  appears  to  give  a 
favourable  comparison  against  the  observed  data.  Reducing  the  bathymetric  data  has  a  major 
influence  on  the  spatial  predictions.  In  some  cases  this  manifests  itself  even  when  half  the 
data  is  used  such  as  at  sites  1,3  and  6,  where  the  spatial  match  degrades  rapidly.  As  the  data 
reduction  continues  huge  errors  may  be  introduced. 


6.3  Summary 

The  results  from  this  section  have  shown  the  bathymetric  data  to  have  a  profound  influence 
on  the  spatial  predictions  of  the  2-D  hydraulic  model  and  a  lesser  influence  on  the  water 
surface  elevation.  This  fits  well  with  the  rigour  of  validation  discussed  in  section  5.  There  it 
was  proposed  that  matching  internal  stage  values  was  a  weaker  form  of  validation  than 
matching  2-D  inundation.  These  results  show  that  a  more  poorly  specified  model,  i.e.  that 
with  less  bathymetric  data,  can  match  the  internal  stage  values  but  not  the  inundation  extent. 
Hence  inundation  extent  validation  is  a  more  rigorous  test  of  model  performance. 

From  these  results  it  would  seem,  as  expected,  that  the  highest  density  of  bathymetric  data 
produces  the  best  results.  More  detailed  simulations  looking  at  the  relationship  between  mesh 
resolution,  bathymetric  data  density  and  distribution,  and  interpolation  method  must  be 
carried  out  before  more  detailed  guidelines  on  the  optimum  level  of  bathymetric  data  can  be 
considered. 


24 


Missouri  Rvuet  Moddlvng  Prcject 


Figure  6-2:  Comparison  of  the  Landsat  TM  image  to  model  predictions  with  the  4  levels  of 
model  bathymetry  at  site  1.  The  bathymetries  were  created  by  varying  the  spacing  of  cross 
sections  used  from  150  m  to  1200  m. 

Figure  6-3:  Comparison  of  the  Landsat  TM  image  to  model  predictions  with  the  4  levels  of 
model  bathymetry  used  at  site  2.  The  bathymetries  were  created  by  varying  the  spacing  of 
cross  sections  used  from  150  m  to  1200  m. 
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Model  Predictions 
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Figure  #=Sa  Comparison  of  the  Landsat  TM  image  to  model  predictions  with  the  4  levels  of 
model  bathymetry  used  for  event  1  at  site  1.  The  bathymetries  were  created  by  varying  the 
spacing  of  cross  sections  used  from  150  m  to  1200  m. 
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Figure  Comparison  of  the  Landsat  TM  image  to  model  predictions  with  the  4  levels  of 
model  bathymetry  used  for  event  1  at  site  2.  The  bathymetries  were  created  by  varying  the 
spacing  of  cross  sections  used  from  150  m  to  1200  m. 
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7.  CONCLUSIONS 


This  modelling  study  of  a  55  km  reach  of  the  Missouri  River  has  illustrated  the  capability  of 
high  resolution,  physically  based  fluvial  hydraulic  models  for  accurate  flow  representation.  A 
hydraulic  model  was  created  for  the  reach  of  interest,  between  Gavins  Point  Dam  and 
Maskell,  using  the  TELEMAC-2D  modelling  system.  A  sensitivity  analysis  of  the  model 
showed  bed  friction  to  be  the  dominant  parameter  governing  model  performance.  Extensive 
validation  of  the  model  predictions  was  carried  out  against  both  internal  stage  data  and 
satellite  imagery  of  flow  field  boundaries.  Finally  the  impact  of  mesh  resolution  and 
bathymetric  data  provision  was  briefly  examined. 

The  key  findings  of  this  study  can  be  summarised  thus: 

•  A  state  of  the  art  two-dimensional  hydraulic  model  can  simulate  fluvial  hydraulics 
accurately  on  this  scale  in  both  steady  state  and  dynamic  conditions. 

•  Model  predictions  can  be  matched  to  stage  data  internal  to  the  model  domain  with  a  high 
degree  of  accuracy  for  initial  simulations  (Figure  5-2).  Improvements  to  the  match  could 
no  doubt  follow  with  further  model  developments. 

•  Spatial  distribution  of  flow  field  boundaries  has  been  shown  to  be  accurate  in  most  cases 
by  comparison  to  satellite  imagery  (Landsat  TM)  (Figure  5-4  and  Figure  5-5). 

•  Partial  model  validation  has  been  graded  by  its  rigour.  Matching  2-D  inundation  data  is 
considered  more  powerful  than  matching  internal  stage  data,  a  hypothesis  corroborated  by 
model  predictions  with  varying  specifications  of  bathymetry  (Figure  6-2,  Figure  6-3  and 
Table  6-1). 

•  A  degree  of  process  validation  was  achieved  against  satellite  imagery  in  a  specific  case 
(Figure  5-6). 

•  The  relationship  between  mesh  resolution  and  bathymetric  data  is  vital  to  model 
specification  but  requires  further  study  to  quantify. 
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8.  RECCOMENDATIONS 

There  remain  however,  many  areas  where  further  work  could  improve  the  model  predictions. 

•  The  fundamental  complexity  and  connection  of  mesh  and  bathymetry  is  a  limiting  factor 
whilst  a)  only  limited  bathymetric  data  is  available  and  b)  computational  expense  must  be 
considered  when  making  meshes.  Until  meshes  can  be  much  more  finely  resolved  and  the 
bathymetric  data  available  to  apply  to  them  at  a  commensurate  scale  then  this  problem  of 
bathymetric  representation  and  its  impact  on  model  predictions  shall  remain.  Further  data 
collection  in  specific  problem  areas  could  mitigate  the  problem  in  this  specific  case. 

•  The  model  set-up  used  here  in  terms  of  parameters  was  extremely  simple  and  as  such 
model  results  were  not  as  good  as  they  could  be  following  further  modifications.  With 
modification  of  the  bed  friction  parameter,  distributing  it  over  time  and/or  space  the  errors 
particularly  in  internal  stage  could  be  reduced  dramatically.  The  2-D  flow  field 
boundaries  are  however  more  dependant  on  the  bathymetry  than  the  parameter 
distribution. 

•  With  the  knowledge  gained  from  the  simulations  carried  out  already  the  model  structure 
could  be  modified  to  include  more  permanent  islands  where  they  are  evident  and  vary  the 
mesh  resolution  over  the  domain  where  the  gradient  of  predicted  variables  are  greatest. 

•  Now  that  model  validation  has  been  somewhat  structured  the  possibility  exists  for 
deliberate  data  collection  in  order  to  validate  specific  outputs  of  the  model  to  produce  as 
strong  and  complete  a  validation  as  possible.  For  example  more  2-D  process  data  could 
be  collected  in  order  to  validate  the  most  fundamental  aspect  of  a  2-D  hydraulic  model. 

The  ever  progressing  capability  of  highly  resolved  2-D  fluvial  hydraulic  models,  including 
flows  over  dry  areas,  opens  up  a  wide  range  of  potential  applications  that  include: 

•  Predicting  the  impact  of  anthropogenic  activity 

•  Flood  hazard  assessment 

•  Geomorphological  and  hydrological  process  studies 

•  River  management  schemes 

This  report  has  shown  that  the  capability  now  exists  for  high  resolution  fluvial  hydraulic 
modelling  to  be  applied  over  long  (>50  km)  reaches  with  a  high  degree  of  accuracy.  The  use 
of  wetting  and  drying  algorithms  allows  the  simulation  of  dynamic  flows  over  initially  dry 
areas  and  their  retreat.  This  capability  opens  up  a  wide  range  of  potential  applications 
especially  for  the  “what  if’  type  modelling  vital  for  controlled  management  and  development 
of  river  systems. 
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Figure  Captions 

Figure  3-J:  The  TELEMAC-2D  modelling  system  (after  Hervouet  and  Lang,  1995).  _ 

Figure  3-2  The  geometry  notation  used  in  TELEMAC-2D  (after  Hervouet  and  Van  Haren,  1995).  _ 

Figure  3-3  Partially  dry  elements:  (a)  the  real  situation,  (b)  element  remaining  fully  wet  introducing  spurious 

momentum  terms,  (c)  element  exclusion  method  and  (d)  TELEMAC-2D  method  keeping  the  element 

partially  dry  (after  Price,  1997). _ _ _ ^ ^ ^  ^ 

Figure  4-la)  Location  and  (b)  detail  of  the  study  reach  of  the  Missouri  River,  USA.  The  modelled  reach  is  from 

Gavins  Point  Dam  to  Masked  and  includes  the  gauging  stations  at  Yankton  and  Gayville. _ 13 

Figure  4-2  The  boundary  delineation  of  the  Missouri  River  model. _ _ 

Figure  4-3  Finite  element  mesh  of  the  Missouri  River  from  Gavins  Point  Dam  to  Masked.  _ ^ - 14 

Figure  4-4  An  example  portion  of  the  models  bathymetric  representation.  This  section  is  immeadiately  upstream 

of  Yankton. _ _ _ _ _ ^ - ; - - — 

Figure  5-1  Linear  sensitivity  of  water  surface  elevation  to  bed  friction  at  the  two  gauge  sites,  Yankton  and 

Gayville. _ _ _ _ _ _ _ _ _ _ _ ^  ^ 

Figure  5-2  Downstream  section  of  the  reach  showing  the  elevations  of  the  river  bed,  observed  water  surface  and 

predicted  water  surface  (using  a  bed  friction  of  n  =  0.025)  at  the  4  gauges.  - ^ - 1^ 

Figure  5-3:  Locations  of  the  two  sites  to  be  used  in  the  comparison  of  model  predictions  with  satellite 

imagery. - - - - ^ - — -  — ^ 

Figure  5-4:  Comparison  ofLandsat  TM  image  and  4  alternative  model  predictions  with  different  bed  friction 

(Manning's  'n')  values  at  site  1.  _ _ _ _ _ _ — - - 

Figure  5-5:  Comparison  ofLandsat  TM  image  and  4  alternative  model  predictions  with  different  bed  friction 

(Manning's  ‘n ')  values  at  site  2.  _ _ _ _ _ ^ 

Figure  5-6  An  overlay  of  the  model  prediction  of flow  field  boundary  onto  the  Landsat  TM  image  for  a  region 
around  Goat  Island  (including  sites  5  and  6)  for  event  1. 


Figure  6-1  Impact  of  degrading  the  number  of  cross  sections  used  on  the  model  bathymetry.  The  spacing 
degrades  from  150  m  to  1200  m  across  the  range  illustrated. 
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Figure  6-2:  Comparison  of  the  Landsat  TM  image  to  model  predictions  with  the  4  levels  of  model  bathymetry 
at  site  1.  The  bathymetries  were  created  by  varying  the  spacing  of  cross  sections  used  from  1 50  m  to  1200  m 

Figure  6-3:  Comparison  of  the  Landsat  TM  image  to  model  predictions  with  the  4  levels  of  model  bathymetry 
used  at  site  2.  The  bathymetries  were  created  by  varying  the  spacing  of  cross  sections  used  from  150  m  to 
1200  m.  _ _ _ _ _ _ _ 


Table  Captions 

Table  4-1  Gauge  station  location  and  data  availability. - - - - 

Table  6-1  Water  surface  elevations  at  Yankton  and  Gayville  for  simulations  with  the  various  levels  of 
bathymetric  data  supplied  to  the  model.  _ _ _ _ _ _ 
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